#include"logdata.h"

static logF* LD_initLog(int data);
static int LD_reallocM(logF *log);
static int LD_reallocS(logF *log);
static int LD_read(FILE* fp,logF *log, int pos[],int nbpos, int nbread, int proc, int bproc);
static int LD_stat(logF* log);

logFinfo* LD_readInfo(int argc,char* argv[]){
  char buffer[PLIB_STR];
  char dummy;
  FILE* logfp;
  char logfile[PLIB_STR];

  logFinfo* info;

  /* alloc */
  info           = U_malloc(1,logFinfo,"LD_readinfo error 1.0\n");
  info->base     = U_malloc(PLIB_STR,char,"LD_readinfo error 1.1\n");
  if(argc == 3)
    info->coords = U_malloc(PLIB_STR,char,"LD_readinfo error 1.2\n");
  else
    info->coords = NULL;
  info->F_graph  = U_malloc(PLIB_STR,char,"LD_readinfo error 1.3\n");
  info->F_grid   = U_malloc(PLIB_STR,char,"LD_readinfo error 1.4\n");

  /* get graph file from log first line */
  strcpy(info->base,argv[1]);
  if(argc == 3)
    strcpy(info->coords,argv[2]);

  sprintf(logfile,"%s.log",argv[1]);
  WARNING(((logfp = fopen(logfile,"r")) == NULL,"LD_readinfo error 2 : cannot open %s\n",logfile),NULL);

  WARNING((fgets(buffer,PLIB_STR,logfp) == NULL,
	   "LD_readinfo error 2 : cannot read line 1\n"),NULL);  
  sscanf(buffer,"%c%s",&dummy,info->F_graph);  
  WARNING((fgets(buffer,PLIB_STR,logfp) == NULL,
	   "LD_readinfo error 3 : cannot read line 2\n"),NULL);  
  sscanf(buffer,"%c%s",&dummy,info->F_grid);

  fclose(logfp);

  return info;
}

void LD_freeInfo(logFinfo* info){
  U_free("LD_freeInfo error 1",&(info->base),&(info->coords),&(info->F_graph),&(info->F_grid),&info,LEND);
}

logF* LD_readTime(FILE *fp, ctrl* ctrl){
  int i;
  logF *log;
  int pos[5],nbpos=5,nbread=8;
  FILE* out;

  switch(ctrl->stat_out){
  case FILES:
    WARNING(((out = fopen(ctrl->stat_file,"w")) == NULL,"LD_readTime error code 1 : cannot open %s\n",ctrl->stat_file),NULL);
    break;
  case STDOUT:
    out = stdout;
    break;
  default:
    NOTICE((1,"LD_readTime error code 1 : unknown stat_out option (%d) take default stdout",ctrl->stat_out));
    out = stdout;
  }
  pos[TIME_COARSENING] = 3;
  pos[TIME_INITIAL] = 4;
  pos[TIME_UNCOARSENING] = 5;
  pos[TIME_EXECTIME] = 6;
  pos[TIME_TOTAL] = 7;
  
  WARNING(((log = LD_initLog(nbpos)) == NULL,"LD_readTime error code 1 : error in LD_initTimeLog\n"),NULL);
  WARNING((LD_read(fp,log,pos,nbpos,nbread,0,0) == 0,"LD_readTime error code 2 : error in LD_read\n"),NULL);
  WARNING((LD_stat(log) == 0,"LD_readTime error code 3 : error in LD_stat\n"),NULL);

  /* output stats */
  fprintf(out,"------------------------------------------------\n------------------------------------------------\nTiming stats\n------------------------------------------------\n");
  for(i=0;i<log->m;i++){
    fprintf(out,"------------------------------------------------\n%s (%d runs)\n",log->methods[i],log->s[i]);
    fprintf(out,"%20s %10s : %20.5lf %20s : %20.5lf\n","coarsening",
	    "average",log->s_data[i][TIME_COARSENING][LD_AVG],
	    "Standard deviation",log->s_data[i][TIME_COARSENING][LD_STDDEV]);
    fprintf(out,"%20s %10s : %20.5lf %20s : %20.5lf\n","initial partitioning",
	    "average",log->s_data[i][TIME_INITIAL][LD_AVG],
	    "Standard deviation",log->s_data[i][TIME_INITIAL][LD_STDDEV]);
    fprintf(out,"%20s %10s : %20.5lf %20s : %20.5lf\n","uncoarsening",
	    "average",log->s_data[i][TIME_UNCOARSENING][LD_AVG],
	    "Standard deviation",log->s_data[i][TIME_UNCOARSENING][LD_STDDEV]);
    fprintf(out,"%20s %10s : %20.5lf %20s : %20.5lf\n","exec time load bal",
	    "average",log->s_data[i][TIME_EXECTIME][LD_AVG],
	    "Standard deviation",log->s_data[i][TIME_EXECTIME][LD_STDDEV]);
    fprintf(out,"%20s %10s : %20.5lf %20s : %20.5lf\n","total",
	    "average",log->s_data[i][TIME_TOTAL][LD_AVG],
	    "Standard deviation",log->s_data[i][TIME_TOTAL][LD_STDDEV]);
  }
  if(out != stdout)
    fclose(out);
  return log;
}

logF* LD_readLog(FILE *fp, ctrl* ctrl){
  int i;
  logF *log;
  int pos[8],nbpos=8,nbread=10;
  FILE* out;
  
  /* options */
  switch(ctrl->stat_out){
  case FILES:
    WARNING(((out = fopen(ctrl->stat_file,"w")) == NULL,"LD_readTime error code 1 : cannot open %s\n",ctrl->stat_file),NULL);
    break;
  case STDOUT:
    out = stdout;
    break;
  default:
    NOTICE((1,"LD_readTime error code 1 : unknown stat_out option take default stdout"));
    out = stdout;
  }

  /* read and do stats */
  pos[LOG_CUTEDGE]              = 3;
  pos[LOG_TOTCOMMCOST]          = 4;
  pos[LOG_MAXWEIGHT]            = -1; /* max first processor series */
  pos[LOG_CUTEDGEBYTOTCOMMCOST] = 5;
  pos[LOG_IMBALANCE]            = 6;
  pos[LOG_MAXEXECTIME]          = 7;
  pos[LOG_TOTEXECTIME]          = 8;
  pos[LOG_IMBEXECTIME]          = 9;  
  
  WARNING(((log = LD_initLog(nbpos)) == NULL,"LD_readLog error code 1 : error in LD_initLog\n"),NULL);
  
  WARNING((LD_read(fp,log,pos,nbpos,nbread,2,5) == 0,"LD_readLog error code 2 : error in LD_read\n"),NULL);
  
  WARNING((LD_stat(log) == 0,"LD_readLog error code 3 : error in LD_stat\n"),NULL);
    
  /* output stats */
  fprintf(out,"------------------------------------------------\n------------------------------------------------\nLog stats\n------------------------------------------------\n");
  for(i=0;i<log->m;i++){
    fprintf(out,"------------------------------------------------\n%s (%d runs)\n",log->methods[i],log->s[i]);
    fprintf(out,"%20s %10s : %20.5lf %20s : %20.5lf\n","cut edge weight",
	    "average",log->s_data[i][LOG_CUTEDGE][LD_AVG],
	    "Standard deviation",log->s_data[i][LOG_CUTEDGE][LD_STDDEV]);
    fprintf(out,"%20s %10s : %20.5lf %20s : %20.5lf\n","",
	    "minumum",log->s_data[i][LOG_CUTEDGE][LD_MIN],
	    "maximum",log->s_data[i][LOG_CUTEDGE][LD_MAX]);
    fprintf(out,"%20s %10s : %20.5lf %20s : %20.5lf\n","total comm cost",
	    "average",log->s_data[i][LOG_TOTCOMMCOST][LD_AVG],
	    "Standard deviation",log->s_data[i][LOG_TOTCOMMCOST][LD_STDDEV]);
    fprintf(out,"%20s %10s : %20.5lf %20s : %20.5lf\n","max weight",
	    "average",log->s_data[i][LOG_MAXWEIGHT][LD_AVG],
	    "Standard deviation",log->s_data[i][LOG_MAXWEIGHT][LD_STDDEV]);
    fprintf(out,"%20s %10s : %20.5lf %20s : %20.5lf\n","CEW / TCC",
	    "average",log->s_data[i][LOG_CUTEDGEBYTOTCOMMCOST][LD_AVG],
	    "Standard deviation",log->s_data[i][LOG_CUTEDGEBYTOTCOMMCOST][LD_STDDEV]);
    fprintf(out,"%20s %10s : %20.5lf %20s : %20.5lf\n","imbalance",
	    "average",log->s_data[i][LOG_IMBALANCE][LD_AVG],
	    "Standard deviation",log->s_data[i][LOG_IMBALANCE][LD_STDDEV]);
    fprintf(out,"%20s %10s : %20.5lf %20s : %20.5lf\n","max exec time",
	    "average",log->s_data[i][LOG_MAXEXECTIME][LD_AVG],
	    "Standard deviation",log->s_data[i][LOG_MAXEXECTIME][LD_STDDEV]);
    fprintf(out,"%20s %10s : %20.5lf %20s : %20.5lf\n","",
	    "minumum",log->s_data[i][LOG_MAXEXECTIME][LD_MIN],
	    "maximum",log->s_data[i][LOG_MAXEXECTIME][LD_MAX]);
    fprintf(out,"%20s %10s : %20.5lf %20s : %20.5lf\n","tot exec time",
	    "average",log->s_data[i][LOG_TOTEXECTIME][LD_AVG],
	    "Standard deviation",log->s_data[i][LOG_TOTEXECTIME][LD_STDDEV]);
    fprintf(out,"%20s %10s : %20.5lf %20s : %20.5lf\n","",
	    "minumum",log->s_data[i][LOG_TOTEXECTIME][LD_MIN],
	    "maximum",log->s_data[i][LOG_TOTEXECTIME][LD_MAX]);
    fprintf(out,"%20s %10s : %20.5lf %20s : %20.5lf\n","imbalance exec time",
	    "average",log->s_data[i][LOG_IMBEXECTIME][LD_AVG],
	    "Standard deviation",log->s_data[i][LOG_IMBEXECTIME][LD_STDDEV]);
    fprintf(out,"%20s %10s : %20.5lf %20s : %20.5lf\n","",
	    "minumum",log->s_data[i][LOG_IMBEXECTIME][LD_MIN],
	    "maximum",log->s_data[i][LOG_IMBEXECTIME][LD_MAX]);
  }

  if(out != stdout)
    fclose(out);

  return log;
}
logF* LD_readBetLog(FILE *fp, ctrl* ctrl){
  int i;
  logF *log;
  int pos[9],nbpos=9,nbread=11;
  FILE* out;
  
  /* options */
  switch(ctrl->stat_out){
  case FILES:
    WARNING(((out = fopen(ctrl->stat_file,"w")) == NULL,"LD_readTime error code 1 : cannot open %s\n",ctrl->stat_file),NULL);
    break;
  case STDOUT:
    out = stdout;
    break;
  default:
    NOTICE((1,"LD_readTime error code 1 : unknown stat_out option take default stdout"));
    out = stdout;
  }
  
  pos[BET_CUTEDGE] = 3;
  pos[BET_TOTCOMMCOST] = 4;
  pos[BET_MAXWEIGHT] = -1; /* max first processor series */
  pos[BET_CUTEDGEBYTOTCOMMCOST] = 5;
  pos[BET_IMBALANCE] = 6;
  pos[BET_MAXEXECTIME_025] = 7;
  pos[BET_MAXEXECTIME_05] = 8;
  pos[BET_MAXEXECTIME_1] = 9;
  pos[BET_MAXEXECTIME_2] = 10;
  
  
  
  WARNING(((log = LD_initLog(nbpos)) == NULL,"LD_readVwtLog error code 1 : error in LD_initTimeLog\n"),NULL);
  
  WARNING((LD_read(fp,log,pos,nbpos,nbread,1,5) == 0,"LD_readVwtLog error code 2 : error in LD_read\n"),NULL);
  
  WARNING((LD_stat(log) == 0,"LD_readVwtLog error code 3 : error in LD_stat\n"),NULL);
    
  /* output stats */
  fprintf(out,"------------------------------------------------\n------------------------------------------------\nVwt Log stats\n------------------------------------------------\n");
  for(i=0;i<log->m;i++){
    fprintf(out,"------------------------------------------------\n%s (%d runs)\n",log->methods[i],log->s[i]);
    fprintf(out,"%20s %10s : %20.5lf %20s : %20.5lf\n","cut edge weight",
	    "average",log->s_data[i][BET_CUTEDGE][LD_AVG],
	    "Standard deviation",log->s_data[i][BET_CUTEDGE][LD_STDDEV]);
    fprintf(out,"%20s %10s : %20.5lf %20s : %20.5lf\n","total comm cost",
	    "average",log->s_data[i][BET_TOTCOMMCOST][LD_AVG],
	    "Standard deviation",log->s_data[i][BET_TOTCOMMCOST][LD_STDDEV]);
    fprintf(out,"%20s %10s : %20.5lf %20s : %20.5lf\n","max weight",
	    "average",log->s_data[i][BET_MAXWEIGHT][LD_AVG],
	    "Standard deviation",log->s_data[i][BET_MAXWEIGHT][1]);
    fprintf(out,"%20s %10s : %20.5lf %20s : %20.5lf\n","CEW / TCC",
	    "average",log->s_data[i][BET_CUTEDGEBYTOTCOMMCOST][LD_AVG],
	    "Standard deviation",log->s_data[i][BET_CUTEDGEBYTOTCOMMCOST][LD_STDDEV]);
    fprintf(out,"%20s %10s : %20.5lf %20s : %20.5lf\n","imbalance",
	    "average",log->s_data[i][BET_IMBALANCE][LD_AVG],
	    "Standard deviation",log->s_data[i][BET_IMBALANCE][LD_STDDEV]);
    fprintf(out,"%20s %10s : %20.5lf %20s : %20.5lf\n","max time RT=0.25",
	    "average",log->s_data[i][BET_MAXEXECTIME_025][LD_AVG],
	    "Standard deviation",log->s_data[i][BET_MAXEXECTIME_025][LD_STDDEV]);
    fprintf(out,"%20s %10s : %20.5lf %20s : %20.5lf\n","max time RT=0.5",
	    "average",log->s_data[i][BET_MAXEXECTIME_05][LD_AVG],
	    "Standard deviation",log->s_data[i][BET_MAXEXECTIME_05][LD_STDDEV]);
    fprintf(out,"%20s %10s : %20.5lf %20s : %20.5lf\n","max time RT=1",
	    "average",log->s_data[i][BET_MAXEXECTIME_1][LD_AVG],
	    "Standard deviation",log->s_data[i][BET_MAXEXECTIME_1][LD_STDDEV]);
    fprintf(out,"%20s %10s : %20.5lf %20s : %20.5lf\n","max time RT=2",
	    "average",log->s_data[i][BET_MAXEXECTIME_2][LD_AVG],
	    "Standard deviation",log->s_data[i][BET_MAXEXECTIME_2][LD_STDDEV]);
  }
  
  if(out != stdout)
    fclose(out);

  return log;
}

logF* LD_initLog(int data){
  logF* log;
  int i;
  
  log = U_malloc(1,logF,"LD_initLog error code 1\n");
  log->d     = data;
  log->max_d = data;
  log->max_m = METHODS;
  log->max_s = SAMPLES;
  log->max_c = COARSE;
  
  /* methods */
  log->m = 0;
  log->methods = U_malloc(METHODS,char *,"LD_initLog error code 2\n");
  
  /* number of samples */
  log->s = U_malloc(METHODS,int,"LD_initLog error code 3\n");
  for(i=0;i<METHODS;i++){
    log->s[i] = 0;
  }

  /* coarsening files */
  log->c = 0;
  log->coarse = U_malloc(log->max_c,char *,"LD_initLog error code 4\n");

  /* values to be readen*/
  log->matching = U_mallocxy(log->max_m,log->max_s,int,"LD_initLog error code 5\n");
  log->id = U_mallocxy(log->max_m,log->max_s,int,"LD_initLog error code 6\n");
  log->data = U_malloc(log->max_m,double **,"LD_initLog error code 7\n");
  log->s_data = U_malloc(log->max_m,double **,"LD_initLog error code 8\n");

  for(i=0;i<log->max_m;i++){
    log->data[i] = U_mallocxy(log->max_d,log->max_s,double,"LD_initLog error code 8\n");
    log->s_data[i] = U_mallocxy(log->max_d,4,double,"LD_initLog error code 9\n");
  }
  /* sums */
  log->sum = U_mallocxy(log->max_m,log->max_d,double,"LD_initLog error code 10\n");
  
  return log;
}

int LD_reallocM(logF *log){
  int old_max_m = log->max_m,i;
  
  log->max_m *= 2;

  log->methods = U_realloc(log->methods,log->max_m,char *,"LD_reallocM error code 1");
  log->sum = U_reallocxy(log->sum,old_max_m,log->d,log->max_m,log->d,double,"LD_reallocM error code 2\n");
  log->matching = U_reallocxy(log->matching,old_max_m,log->max_s,log->max_m,log->max_s,int,"LD_reallocM error code 3\n");
  log->id = U_reallocxy(log->id,old_max_m,log->max_s,log->max_m,log->max_s,int,"LD_reallocM error code 3\n");
  log->data = U_realloc (log->data,log->max_m,double **,"LD_reallocM error code 4\n");
  log->s_data = U_realloc (log->s_data,log->max_m,double **,"LD_reallocM error code 5\n");
  for(i=old_max_m;i<log->max_m;i++){
    log->data[i] = U_mallocxy(log->max_d,log->max_s,double,"LD_reallocM error code 6\n");
    log->s_data[i] = U_mallocxy(log->max_d,4,double,"LD_reallocM error code 7\n");
  }
  log->s = U_realloc(log->s,log->max_m,int,"LD_reallocM error code 7");
  for(i=old_max_m;i<log->max_m;i++)
    log->s[i] = 0;

  return 1;
}

int LD_reallocS(logF *log){
  int i;
  int 
    x     = log->max_m,
    y     = 2*log->max_s,
    old_x = log->max_m,
    old_y = log->max_s;
  printf("reallocS\n");
  log->max_s*=2;
  log->matching = U_reallocxy(log->matching,old_x,old_y,x,y,int,
			      "LD_reallocS error code 1\n");
  log->id       = U_reallocxy(log->id,old_x,old_y,x,y,int,
			      "LD_reallocS error code 2\n");
  x = old_x = log->max_d;
  for(i=0;i<log->max_m;i++)
    log->data[i] = U_reallocxy(log->data[i],old_x,old_y,x,y,double,
			       "LD_reallocS error code 3\n");
  return 1;
}

int LD_read(FILE* fp,logF *log, int pos[],int nbpos, int nbread, int proc, int bproc){
  int read,l=0,nbproc = 0;
  double max = 0.0;
  int i,j,k;
  
  while(!feof(fp)){
    /* read line */
    l++;
    /* check that 3rd field is not alpha or empty in such a case ignore line */    
    if(((read = U_readLine(fp)) == 0) || (isalpha(U_line[3][0])) || (strlen(U_line[3])==0)){
      /* if header line count the number of processors */
      if(isalpha(U_line[3][0])){
	if(proc != 0)
	  nbproc = (read-nbread)/proc;
	else
	  nbproc = 0;
      }
      continue;
    }
    WARNING((read != nbread + proc * nbproc,"LD_read error code 1 : error line %d in log file (read %d, expected %d (%d processors)\n",l,read,nbread + proc * nbproc,nbproc),0);
    
    /* get method id (i) */
    for(i=0;i<log->m;i++){
      if(strcmp(log->methods[i],U_line[2]) == 0)
	break;
    }
    if(i == log->m){
      log->methods[i] = U_malloc(strlen(U_line[2])+1,char,"LD_read error code 2\n");
      strcpy(log->methods[i],U_line[2]);
      log->m++;
      if(log->m >= log->max_m)
	WARNING((LD_reallocM(log) == 0,"LD_read error code 3 : error in LD_reallocM\n"),0);
    }
    /* put matching serialize file into matching */
    if(strlen(U_line[1]) == 0)
      j = -1;
    else{
      for(j=0;j<log->c;j++){
	if(strcmp(log->coarse[j],U_line[1]) == 0)
	  break;
      }
    }
    if(j == log->c){
      log->coarse[j] = U_malloc(strlen(U_line[1])+1,char,"LD_read error code 4\n");
      strcpy(log->coarse[j],U_line[1]);
      log->c++;
      if(log->c >= log->max_c){
	log->max_c *= 2;
	log->coarse = U_realloc(log->coarse,log->max_c,char *,"LD_read error code 5\n");
      }
    }

    log->matching[i][log->s[i]] = j;

    /* result id */
    log->id[i][log->s[i]] = atoi(U_line[0]);
    
    /* put fields into data */
    for(j=0;j<nbpos;j++){
      max = 0.0;
      if(pos[j]<0){
	for(k = bproc + 1 + ((-pos[j]-1) * nbproc);k<bproc + 1 + ((-pos[j]-1) * nbproc) + nbproc ; k++){
	  if(atof(U_line[k]) > max)
	    max = atof(U_line[k]);
	}
	log->data[i][j][log->s[i]] = max;
      }else if(pos[j]<=bproc){
	log->data[i][j][log->s[i]] = atof(U_line[pos[j]]);
      }else{
	/*printf("[ %d %d %d ] %d %s %d %s %d %s\n",pos[j],proc,nbproc,
	       pos[j] + proc * nbproc-2,U_line[pos[j] + proc * nbproc-2],
	       pos[j] + proc * nbproc-1,U_line[pos[j] + proc * nbproc-1],
	       pos[j] + proc * nbproc,U_line[pos[j] + proc * nbproc]);*/
	log->data[i][j][log->s[i]] = atof(U_line[pos[j] + proc * nbproc]);
      }
      
    }
    
    log->s[i]++;
    if(log->s[i] >= log->max_s)
      WARNING((LD_reallocS(log) == 0,"LD_read error code 6 : error in LD_reallocS\n"),0);
  }  
  
  U_freeline();
  return 1;
}

int LD_stat(logF* log){
  int i,j,k;
  double min,max;

  for(i=0;i<log->m;i++){
    for(k=0;k<log->d;k++){
      log->sum[i][k] = 0.0;

      min = +INFINITY;
      max = -INFINITY;
      for(j=0;j<log->s[i];j++){

	if(log->data[i][k][j] > max)
	  max = log->data[i][k][j];
	if(log->data[i][k][j] < min)
	  min = log->data[i][k][j];

	log->sum[i][k] += log->data[i][k][j];
      }
      
      log->s_data[i][k][LD_AVG] = log->sum[i][k] / log->s[i];
      log->s_data[i][k][LD_MIN] = min;
      log->s_data[i][k][LD_MAX] = max;
    }
  }
  for(i=0;i<log->m;i++){
    for(k=0;k<log->d;k++){
      log->sum[i][k] = 0.0;
      for(j=0;j<log->s[i];j++)
	log->sum[i][k] += SQUARE(log->data[i][k][j] - log->s_data[i][k][0]);
      log->s_data[i][k][LD_STDDEV] = sqrt(log->sum[i][k] / log->s[i]);
    }
  }
  return 1;
}

double** LD_getData(logF *log,int x, int *methods,int nbmethods,int group, int* nbdat){
  int i, j, k, l;
  double **data;
  int s;

  if(group == BY_COARSENING){
    if(methods == NULL || nbmethods == 0){
      *nbdat  = log->m;
      data = U_malloc(log->m,double *,"LD_getData error code 1\n");
      for(i=0;i<log->m;i++){
	/* get matching data */
	data[i] = U_malloc(log->s[i]+1,double,"LD_getData error code 2\n");
	data[i][0] = log->s[i];
	memcpy(&(data[i][1]),log->data[i][x],log->s[i]*sizeof(double));
      }
    }else{
      *nbdat  = nbmethods;
      data = U_malloc(nbmethods,double *,"LD_getData error code 3\n");
      for(j=0;j<nbmethods;j++){
	i = methods[j];
	/* get matching data */
	data[j] = U_malloc(log->s[i]+1,double,"LD_getData error code 4\n");
	data[j][0] = log->s[i];
	memcpy(&(data[j][1]),log->data[i][x],log->s[i]*sizeof(double));
      }
    }
    return data;
  }
  
  if(group == BY_MATCHING){
    if(methods == NULL || nbmethods == 0){
      *nbdat = log->c;
      data = U_mallocxy(log->c,log->max_s+1,double,"LD_getData error code 5\n");
      for(i=0;i<log->c;i++){
	s = 0;
	for(j=0;j<log->m;j++){
	  for(k=0;k<log->s[j];k++){
	    if(log->matching[j][k] == i){
	      data[i][s+1] = log->data[j][x][k];
	      s++;
	    } 
	  }
	}
	data[i][0] = s;
      }
    }else{
      *nbdat = nbmethods;
      data = U_mallocxy(nbmethods,log->max_s+1,double,"LD_getData error code 7\n");
      for(l=0;l<nbmethods;l++){
	i = methods[l];
	s = 0;
	for(j=0;j<log->m;j++){
	  for(k=0;k<log->s[j];k++){
	    if(log->matching[j][k] == i){
	      data[l][s+1] = log->data[j][x][k];
	      s++;
	    } 
	  }
	}
	data[l][0] = s;
      }
    }
    return data;
  }
  return NULL;
}

plot *LD_getPlotXY(logF *log,int x,int y,int *methods,int nbmethods,int group,char *xlabel,char *ylabel, char *title, int border, int grid, int key, int xtics, int ytics){
  plot *p;
  data *dat;
  int i,j,s=0,k,l;

  p = U_malloc (1,plot,"LD_getPlotXY error code 1\n");

  p->xlabel = xlabel;
  p->ylabel = ylabel;
  p->title  = title;
  p->border = border;
  p->grid   = grid;
  p->key    = key;
  p->xtics  = xtics;
  p->ytics  = ytics;
  p->wins   = NULL;

  if(group == BY_COARSENING){
    dat = U_malloc(log->m,data,"LD_getPlotXY error code 2\n");
    p->nbdat  = nbmethods;
    for(j=0;j<nbmethods;j++){
      i = methods[j];
      dat[j].pdata = NULL;
      dat[j].nb = log->s[i];
      dat[j].x  = log->data[i][x];
      dat[j].y  = log->data[i][y];
      dat[j].z  = NULL;
      dat[j].title = log->methods[i];
      dat[j].style = PLOT_POINTS;
    }
    p->dat=dat;
    return p;
  }

  if(group == BY_MATCHING){
    dat = U_malloc(log->c,data,"LD_getPlotXY error code 3\n");
    for(l=0;l<nbmethods;l++){
      i = methods[l];
      dat[l].pdata = U_malloc(2 * log->max_s,double,"LD_getPlotXY error code 4\n");
      dat[l].x = dat[l].pdata;
      dat[l].y = dat[l].pdata + log->max_s;
      dat[l].z = NULL;
      s = 0;
      for(j=0;j<log->m;j++){
	for(k=0;k<log->s[j];k++){
	  if(log->matching[j][k] == i){
	    dat[l].x[s] = log->data[j][x][k];
	    dat[l].y[s] = log->data[j][y][k];
	    s++;
	  } 
	}
      }
      dat[l].nb = s;
      dat[l].title = log->coarse[i];
      dat[l].style = PLOT_POINTS;
    }
    p->nbdat = nbmethods;
    p->dat = dat;
    return p;
  }
  return NULL;
}
plot *LD_getPlotXYZ(logF *log,int x,int y,int z,int *methods,int nbmethods,int group,char *xlabel,char *ylabel, char* zlabel, char *title, int border, int grid, int key, int xtics, int ytics){
  plot *p;
  data *dat;
  int i,j,s=0,k,l;

  p = U_malloc (1,plot,"LD_getPlotXY error code 1\n");

  p->xlabel = xlabel;
  p->ylabel = ylabel;
  p->zlabel = zlabel;
  p->title  = title;
  p->border = border;
  p->grid   = grid;
  p->key    = key;
  p->xtics  = xtics;
  p->ytics  = ytics;
  p->wins   = NULL;

  if(group == BY_COARSENING){
    dat = U_malloc(log->m,data,"LD_getPlotXY error code 2\n");
    p->nbdat  = nbmethods;
    for(j=0;j<nbmethods;j++){
      i = methods[j];
      dat[j].pdata = NULL;
      dat[j].nb = log->s[i];
      dat[j].x  = log->data[i][x];
      dat[j].y  = log->data[i][y];
      dat[j].z  = log->data[i][z];
      dat[j].title = log->methods[i];
      dat[j].style = PLOT_POINTS;
    }
    p->dat=dat;
    return p;
  }

  if(group == BY_MATCHING){
    dat = U_malloc(log->c,data,"LD_getPlotXY error code 3\n");
    for(l=0;l<nbmethods;l++){
      i = methods[l];
      dat[l].pdata = U_malloc(3 * log->max_s,double,"LD_getPlotXY error code 4\n");
      dat[l].x = dat[l].pdata;
      dat[l].y = dat[l].pdata + log->max_s;
      dat[l].z = dat[l].pdata + 2*log->max_s;
      s = 0;
      for(j=0;j<log->m;j++){
	for(k=0;k<log->s[j];k++){
	  if(log->matching[j][k] == i){
	    dat[l].x[s] = log->data[j][x][k];
	    dat[l].y[s] = log->data[j][y][k];
	    dat[l].z[s] = log->data[j][z][k];
	    s++;
	  } 
	}
      }
      dat[l].nb = s;
      dat[l].title = log->coarse[i];
      dat[l].style = PLOT_POINTS;
    }
    p->nbdat = nbmethods;
    p->dat = dat;
    return p;
  }
  return NULL;
}

double** LD_getDataMatch(logF *log,double *x, int *methods,int nbmethods,int group, int* nbdat){
  int i, j, k;
  double **data;

  if(group == BY_COARSENING){
    if(methods == NULL || nbmethods == 0){
      *nbdat  = log->m;
      data = U_malloc(log->m,double *,"LD_getDataMatch error code 1\n");
      for(i=0;i<log->m;i++){
	/* get matching data */
	data[i] = U_malloc(log->s[i]+1,double,"LD_getDataMatch error code 2\n");
	data[i][0] = log->s[i];
	for(k=0;k<log->s[i];k++){
	  data[i][k+1] = x[log->matching[i][k]];
	}
      }
    }else{
      *nbdat  = nbmethods;
      data = U_malloc(nbmethods,double *,"LD_getDataMatch error code 3\n");
      for(j=0;j<nbmethods;j++){
	i = methods[j];
	/* get matching data */
	data[j] = U_malloc(log->s[i]+1,double,"LD_getPlotXYmatchX error code 3\n");
	data[j][0] = log->s[i];
	for(k=0;k<log->s[i];k++){
	  data[j][k+1] = x[log->matching[i][k]];
	}
      }
    }
    return data;
  }
  return NULL;
}

plot *LD_getPlotXYmatchX(logF *log,int y,double* x,int *methods,int nbmethods,int group,char *xlabel,char *ylabel, char *title, int border, int grid, int key, int xtics, int ytics){
  plot *p;
  data *dat;
  int i,j,k,s=0,l;
  double *m_x;

  p = U_malloc (1,plot,"LD_getPlot error code 1\n");

  p->xlabel = xlabel;
  p->ylabel = ylabel;
  p->title  = title;
  p->border = border;
  p->grid   = grid;
  p->key    = key;
  p->xtics  = xtics;
  p->ytics  = ytics;
  p->wins   = NULL;

  if(group == BY_COARSENING){
    dat = U_malloc (log->m,data,"LD_getPlot error code 2 : memory allocation error\n");
    p->nbdat  = 0;
    for(j=0;j<nbmethods;j++){
      i = methods[j];
      /* get matching data */
      dat[j].pdata = m_x = U_malloc(log->s[i],double,"LD_getPlotXYmatchX error code 3\n");
      for(k=0;k<log->s[i];k++){
	m_x[k] = x[log->matching[i][k]];
	/*printf("points %lf %lf\n",log->data[i][y][k],m_x[k]);*/
      }
      /* other data (from log files) */
      dat[j].nb = log->s[i];
      dat[j].x  = m_x;
      dat[j].y  = log->data[i][y];
      dat[j].z  = NULL;
      dat[j].title = log->methods[i];
      dat[j].style = PLOT_POINTS;
    }
    p->nbdat = nbmethods;
    p->dat=dat;
    return p;
  }
  if(group == BY_MATCHING){
    dat = U_malloc(log->c,data,"LD_getPlotXY error code 3\n");
    for(l=0;l<nbmethods;l++){
      i = methods[l];
      dat[l].pdata = U_malloc(2 * log->max_s,double,"LD_getPlotXY error code 4\n");
      dat[l].x = dat[l].pdata;
      dat[l].y = dat[l].pdata + log->max_s;
      dat[l].z = NULL;
      s = 0;
      for(j=0;j<log->m;j++){
	for(k=0;k<log->s[j];k++){
	  if(log->matching[j][k] == i){
	    dat[l].x[s] = x[log->matching[j][k]];
	    dat[l].y[s] = log->data[j][y][k];
	    s++;
	  } 
	}
      }
      dat[l].nb = s;
      dat[l].title = log->coarse[i];
      dat[l].style = PLOT_POINTS;
    }
    p->nbdat = nbmethods;
    p->dat = dat;
  }
  return p;
}

plot *LD_getPlotXYZmatchX(logF *log,int y,int z,double* x,int *methods,int nbmethods,int group,char *xlabel,char *ylabel,char *zlabel,char *title, int border, int grid, int key, int xtics, int ytics){
  plot *p;
  data *dat;
  int i,j,k,s=0,l;
  double *m_x;

  p = U_malloc (1,plot,"LD_getPlot error code 1\n");

  p->xlabel = xlabel;
  p->ylabel = ylabel;
  p->zlabel = zlabel;
  p->title  = title;
  p->border = border;
  p->grid   = grid;
  p->key    = key;
  p->xtics  = xtics;
  p->ytics  = ytics;
  p->wins   = NULL;

  if(group == BY_COARSENING){
    dat = U_malloc (log->m,data,"LD_getPlot error code 2 : memory allocation error\n");
    p->nbdat = nbmethods;
    for(j=0;j<nbmethods;j++){
      i = methods[j];
      /* get matching data */
      dat[j].pdata = m_x = U_malloc(log->s[i],double,"LD_getPlotXYmatchX error code 3\n");
      for(k=0;k<log->s[i];k++){
	m_x[k] = x[log->matching[i][k]];
	/*printf("points %lf %lf\n",log->data[i][y][k],m_x[k]);*/
      }
      /* other data (from log files) */
      dat[j].nb = log->s[i];
      dat[j].x  = m_x;
      dat[j].y  = log->data[i][y];
      dat[j].z  = log->data[i][z];
      dat[j].title = log->methods[i];
      dat[j].style = PLOT_POINTS;
    }   
    p->dat=dat;
    return p;
  }
  if(group == BY_MATCHING){
    dat = U_malloc(log->c,data,"LD_getPlotXY error code 3\n");
    for(l=0;l<nbmethods;l++){
      i = methods[l];
      dat[l].pdata = U_malloc(3 * log->max_s,double,"LD_getPlotXY error code 4\n");
      dat[l].x = dat[l].pdata;
      dat[l].y = dat[l].pdata + log->max_s;
      dat[l].z = dat[l].pdata + 2*log->max_s;
      s = 0;
      for(j=0;j<log->m;j++){
	for(k=0;k<log->s[j];k++){
	  if(log->matching[j][k] == i){
	    dat[l].x[s] = x[log->matching[j][k]];
	    dat[l].y[s] = log->data[j][y][k];
	    dat[l].z[s] = log->data[j][z][k];
	    s++;
	  } 
	}
      }
      dat[l].nb = s;
      dat[l].title = log->coarse[i];
      dat[l].style = PLOT_POINTS;
    }
    p->nbdat = nbmethods;
    p->dat = dat;
  }
  return p;
}
plot *LD_getPlotXYZmatchXY(logF *log,int z,double* x,double* y,int *methods,int nbmethods,int group,char *xlabel,char *ylabel,char *zlabel,char *title, int border, int grid, int key, int xtics, int ytics){
  plot *p;
  data *dat;
  int i,j,k,s=0,l;
  double *m_x;
  double *m_y;

  p = U_malloc (1,plot,"LD_getPlot error code 1\n");

  p->xlabel = xlabel;
  p->ylabel = ylabel;
  p->zlabel = zlabel;
  p->title  = title;
  p->border = border;
  p->grid   = grid;
  p->key    = key;
  p->xtics  = xtics;
  p->ytics  = ytics;
  p->wins   = NULL;

  if(group == BY_COARSENING){
    dat = U_malloc (log->m,data,"LD_getPlot error code 2 : memory allocation error\n");
    p->nbdat = nbmethods;
    for(j=0;j<nbmethods;j++){
      i = methods[j];
      /* get matching data */
      dat[j].pdata = U_malloc(2 * log->s[i],double,"LD_getPlotXYmatchX error code 3\n");
      m_x = dat[j].pdata;
      m_y = m_x + log->s[i];
      for(k=0;k<log->s[i];k++){
	m_x[k] = x[log->matching[i][k]];
	m_y[k] = y[log->matching[i][k]];
	/*printf("points %lf %lf\n",log->data[i][y][k],m_x[k]);*/
      }
      /* other data (from log files) */
      dat[j].nb = log->s[i];
      dat[j].x  = m_x;
      dat[j].y  = m_y;
      dat[j].z  = log->data[i][z];
      dat[j].title = log->methods[i];
      dat[j].style = PLOT_POINTS;
    }   
    p->dat=dat;
    return p;
  }
  if(group == BY_MATCHING){
    dat = U_malloc(log->c,data,"LD_getPlotXY error code 3\n");
    for(l=0;l<nbmethods;l++){
      i = methods[l];
      dat[l].pdata = U_malloc(3 * log->max_s,double,"LD_getPlotXY error code 4\n");
      dat[l].x = dat[l].pdata;
      dat[l].y = dat[l].pdata + log->max_s;
      dat[l].z = dat[l].pdata + 2*log->max_s;
      s = 0;
      for(j=0;j<log->m;j++){
	for(k=0;k<log->s[j];k++){
	  if(log->matching[j][k] == i){
	    dat[l].x[s] = x[log->matching[j][k]];
	    dat[l].y[s] = y[log->matching[j][k]];
	    dat[l].z[s] = log->data[j][z][k];
	    s++;
	  } 
	}
      }
      dat[l].nb = s;
      dat[l].title = log->coarse[i];
      dat[l].style = PLOT_POINTS;
    }
    p->nbdat = nbmethods;
    p->dat = dat;
  }
  return p;
}

plot *LD_getPlotXYiniX(logF *log,int y,double** x,int *methods,int nbmethods,int group,char *xlabel,char *ylabel, char *title, int border, int grid, int key, int xtics, int ytics){
  plot *p;
  data *dat;
  int i,j,k,s=0,l;
  p = U_malloc (1,plot,"LD_getPlot error code 1\n");
  p->xlabel = xlabel;
  p->ylabel = ylabel;
  p->title  = title;
  p->border = border;
  p->grid   = grid;
  p->key    = key;
  p->xtics  = xtics;
  p->ytics  = ytics;
  p->wins   = NULL;

  if(group == BY_COARSENING){
    dat = U_malloc (log->m,data,"LD_getPlot error code 2 : memory allocation error\n");
    for(j=0;j<nbmethods;j++){
      i = methods[j];
      /* get matching data */
      dat[j].pdata = NULL;
      dat[j].nb = log->s[i];
      dat[j].x  = x[j];
      dat[j].y  = log->data[i][y];
      dat[j].z  = NULL;
      dat[j].title = log->methods[i];
      dat[j].style = PLOT_POINTS;
      p->nbdat++;
    }
    p->nbdat  = nbmethods;
    p->dat=dat;
    return p;
  }
  if(group == BY_MATCHING){
    dat = U_malloc(log->c,data,"LD_getPlotXY error code 3\n");
    for(l=0;l<nbmethods;l++){
      i = methods[l];
      dat[l].pdata = U_malloc(log->max_s,double,"LD_getPlotXY error code 4\n");
      dat[l].y = dat[l].pdata;
      s = 0;
      for(j=0;j<log->m;j++){
	for(k=0;k<log->s[j];k++){
	  if(log->matching[j][k] == i){
	    dat[l].y[s] = log->data[j][y][k];
	    s++;
	  } 
	}
      }
      dat[l].nb = s;
      dat[l].x = x[l];
      dat[l].z = NULL;
      dat[l].title = log->coarse[i];
      dat[l].style = PLOT_POINTS;
    }
    p->nbdat = nbmethods;
    p->dat = dat;
  }
  return p;
}

plot *LD_getPlotXYZiniX(logF *log,int y,int z,double** x,int *methods,int nbmethods,int group,char *xlabel,char *ylabel,char *zlabel,char *title, int border, int grid, int key, int xtics, int ytics){
  plot *p;
  data *dat;
  int i,j,k,s=0,l;
  p = U_malloc (1,plot,"LD_getPlot error code 1\n");
  p->xlabel = xlabel;
  p->ylabel = ylabel;
  p->zlabel = zlabel;
  p->title  = title;
  p->border = border;
  p->grid   = grid;
  p->key    = key;
  p->xtics  = xtics;
  p->ytics  = ytics;
  p->wins   = NULL;

  if(group == BY_COARSENING){
    dat = U_malloc (log->m,data,"LD_getPlot error code 2 : memory allocation error\n");
    for(j=0;j<nbmethods;j++){
      i = methods[j];
      /* get matching data */
      dat[j].pdata = NULL;
      dat[j].nb = log->s[i];
      dat[j].x  = x[j];
      dat[j].y  = log->data[i][y];
      dat[j].z  = log->data[i][z];
      dat[j].title = log->methods[i];
      dat[j].style = PLOT_POINTS;
      p->nbdat++;
    }
    p->nbdat  = nbmethods;
    p->dat=dat;
    return p;
  }
  if(group == BY_MATCHING){
    dat = U_malloc(log->c,data,"LD_getPlotXY error code 3\n");
    for(l=0;l<nbmethods;l++){
      i = methods[l];
      dat[l].pdata = U_malloc(2*log->max_s,double,"LD_getPlotXY error code 4\n");
      dat[l].y = dat[l].pdata;
      dat[l].z = dat[l].pdata + log->max_s;
      s = 0;
      for(j=0;j<log->m;j++){
	for(k=0;k<log->s[j];k++){
	  if(log->matching[j][k] == i){
	    dat[l].y[s] = log->data[j][y][k];
	    dat[l].z[s] = log->data[j][z][k];
	    s++;
	  } 
	}
      }
      dat[l].nb = s;
      dat[l].x = x[l];
      dat[l].title = log->coarse[i];
      dat[l].style = PLOT_POINTS;
    }
    p->nbdat = nbmethods;
    p->dat = dat;
  }
  return p;
}

plot *LD_getPlotXYZiniXY(logF *log,int z,double** x,double** y,int *methods,int nbmethods,int group,char *xlabel,char *ylabel,char *zlabel,char *title, int border, int grid, int key, int xtics, int ytics){
  plot *p;
  data *dat;
  int i,j,k,s=0,l;
  double *ids=NULL;
  double *li[3];

  p = U_malloc (1,plot,"LD_getPlot error code 1\n");
  p->xlabel = xlabel;
  p->ylabel = ylabel;
  p->zlabel = zlabel;
  p->title  = title;
  p->border = border;
  p->grid   = grid;
  p->key    = key;
  p->xtics  = xtics;
  p->ytics  = ytics;
  p->wins   = NULL;

  if(group == BY_COARSENING){
    dat = U_malloc (log->m,data,"LD_getPlot error code 2 : memory allocation error\n");
    for(j=0;j<nbmethods;j++){
      i = methods[j];
      /* get matching data */
      dat[j].pdata = NULL;
      dat[j].nb = log->s[i];
      dat[j].x  = x[j];
      dat[j].y  = y[j];
      dat[j].z  = log->data[i][z];
      for(k=0;k<log->s[i];k++)
	ids[k] = (double)log->id[i][k];
      dat[j].title = log->methods[i];
      dat[j].style = PLOT_POINTS;
      p->nbdat++;
      li[0] = dat[j].y;
      li[1] = dat[j].z;
      li[2] = ids;
      //U_quickSort(dat[j].x,0,dat[j].nb-1,li,3);
    }
    p->nbdat  = nbmethods;
    p->dat=dat;
    return p;
  }
  if(group == BY_MATCHING){
    dat = U_malloc(log->c,data,"LD_getPlotXY error code 3\n");
    for(l=0;l<nbmethods;l++){
      i = methods[l];
      dat[l].pdata = U_malloc(2*log->max_s,double,"LD_getPlotXY error code 4\n");
      dat[l].z = dat[l].pdata;
      ids = dat[l].pdata + log->max_s;
      s = 0;
      for(j=0;j<log->m;j++){
	for(k=0;k<log->s[j];k++){
	  if(log->matching[j][k] == i){
	    dat[l].z[s] = log->data[j][z][k];
	    ids[s] = (double)log->id[j][k];
	    s++;
	  } 
	}
      }
      dat[l].nb = s;
      dat[l].x = x[l];
      dat[l].y = y[l];
      dat[l].title = log->coarse[i];
      dat[l].style = PLOT_POINTS;
      li[0] = dat[l].y;
      li[1] = dat[l].z;
      li[2] = ids;
      U_quickSort(dat[l].x,0,dat[l].nb-1,li,3);
      printf("%14s %10s %10s %10s\n","id",xlabel,ylabel,zlabel);
      for(j=0;j<dat[l].nb;j++){
	printf("%3d %5.0f %10.5f %10.2f %10.2f\n",j,ids[j],dat[l].x[j],dat[l].y[j],dat[l].z[j]);
      }
    }
    p->nbdat = nbmethods;
    p->dat = dat;
  }
  return p;
}

void LD_free(logF *log){
  int i;

  if(log != NULL){
    for(i=0;i<log->max_m;i++)
      U_free("LD_free error code 1\n",&(log->data[i]),&(log->s_data[i]),LEND);
      
    U_free("LD_free error code 1\n",&(log->sum),&(log->matching),&(log->id),LEND);
    
    for(i=0;i<log->m;i++)
      U_free("LD_free error code 4\n",&(log->methods[i]),LEND);
    
    for(i=0;i<log->c;i++)
      U_free("LD_free error code 5\n",&(log->coarse[i]),LEND);
    
    U_free("LD_free error code 6\n",&(log->coarse),&(log->data),&(log->s_data),&(log->methods),&(log->s),&(log),LEND);
  }else{
    NOTICE((1,"LD_free error code 14 : attempting to free NULL log data structure\n"));
  }
}
/* U_free(match,val,id) */
/* last result overall best */
int LD_logBestByMethod(logF *l,int f,int s,char*** match,double **val,int** id){
  int b_s = -1;
  int i,j;
  double b=-1, ***data = l->data;

  *match  = U_malloc(l->m + 1,char*,"LD_logBestByMethod error code 2\n");
  *val    = U_malloc(l->m + 1,double,"LD_logBestByMethod error code 3\n");
  if(id != NULL)
    *id     = U_malloc(l->m + 1,int,"LD_logBestByMethod error code 4\n");
  for(i=0;i<l->m;i++){
    b = -1;
    b_s = -1;
    for(j=0;j<l->s[i];j++){
      switch(s){
      case LD_MIN:
	if(data[i][f][j] < b || b == -1){
	  b = data[i][f][j];
	  b_s = j;
	}
	break;
      case LD_MAX:
	if(data[i][f][j] > b){
	  b = data[i][f][j];
	  b_s = j;
	}
	break;
      default:
	WARNING((1,
		 "LD_logBestMatching error code 1 : unknown search method (LD_MIN or LD_MAX)\n"),
		0);
      }
    }
    if(id != NULL)
      (*id)[i] = l->id[i][b_s];
    (*val)[i]   = b;
    if(l->matching[i][b_s] == -1)
      (*match)[i] = NULL;
    else
      (*match)[i] = l->coarse[l->matching[i][b_s]];
  }
  /* search overall best */
  b = -1;
  b_s = -1;
  for(i=0;i<l->m;i++){
    switch(s){
    case LD_MIN:
      if((*val)[i] < b || b == -1){
	b =  (*val)[i];
	b_s = i;
      }
      break;
    case LD_MAX:
      if((*val)[i] > b){
	b =  (*val)[i];
	b_s = i;
      }
      break;
    }
  }
  if(id != NULL)
    (*id)[l->m] = (*id)[b_s];
  (*val)[l->m] = (*val)[b_s];
  (*match)[l->m] = (*match)[b_s];
  return 1;
}
char* LD_logBestMatching(logF *l,int f,int s,double *val,int* id){
  int b_m = -1;
  int b_s = -1;
  double b=-1, ***data = l->data;
  int i,j;
  /* search best matching among each method and samples 
     among field define by m */
  
  for(i=0;i<l->m;i++){
    for(j=0;j<l->s[i];j++){
      switch(s){
      case LD_MIN:
	if(data[i][f][j] < b || b == -1){
	  b = data[i][f][j];
	  b_m = i;
	  b_s = j;
	}
	break;
      case LD_MAX:
	if(data[i][f][j] > b){
	  b = data[i][f][j];
	  b_m = i;
	  b_s = j;
	}
	break;
      default:
	WARNING((1,"LD_logBestMatching error code 1 : unknown search method (LD_MIN or LD_MAX)\n"),NULL);
      }
    }
  }
  *val = b;
  return l->coarse[l->matching[b_m][b_s]];
}

int* LD_getMatchingIds(logF* log,int coarse){
  int 
    max_ids = IDS,
    *ids = U_malloc(IDS,int,"LD_getMatchingIds error code 1\n"),
    nb_ids = 1,
    i,j;

  for(i=0;i<log->m;i++){
    for(j=0;j<log->s[i];j++){
      if(log->matching[i][j] == coarse){
	ids[nb_ids++] = log->id[i][j];
	if(nb_ids>=max_ids){
	  max_ids *= 2;
	  ids = U_realloc(ids,max_ids,int,"LD_getMatchingIds error code 1\n");
	}
      }
    }
  }
  ids[0] = nb_ids - 1;  
  return ids;
}
